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Abstract. We present new results on instabilities in rapidly and differentially rotating 
neutron stars. We model the stars in full general relativity and describe the stellar matter 
adopting a cold realistic equation of state based on the unified SLy prescription (T]. We 
provide evidence that rapidly and differentially rotating stars that are below the expected 
threshold for the dynamical bar-mode instability, /3 C = T/|iy| ~ 0.25, do nevertheless 
develop a shear instability on a dynamical timescale and for a wide range of values of 0. 
This class of instability, which has so far been found only for small values of ft and with 
very small growth rates, is therefore more generic than previously found and potentially more 
effective in producing strong sources of gravitational waves. Overall, our findings support the 
phenomenological predictions made by Watts, Andersson and Jones (2j on the nature of the 
low-T/|Vy | instability as the manifestation of a shear instability in a region where the latter is 
possible only for small values of the fi. Furthermore, our results provide additional insight on 
shear instabilities and on the necessary conditions for their development. 



PACS numbers: 04.40.Dg, 95.30.Lz, 95.30.Sf 97.60.Jd 

1. Introduction 

Non-axisymmetric deformations of rapidly rotating bodies are rather generic phenomena 
in Nature and can appear in a wide class of systems. Particularly interesting within an 
astrophysical context are those deformations taking place in fluids that are self-gravitating 
and the literature on this has a long history dating back to the work of [3) on incompressible 
Newtonian uniformly rotating bodies. Since then, the study of these instabilities has continued 
over the years both in Newtonian gravity and in full general relativity. 

Special attention has traditionally been paid to the study of m = 2 instabilities, which 
are characterized by the exponential growth of m = 2 deformations, where m parametrizes 
the azimuthal dependence e lm ^ in a standard mode decomposition in spherical harmonics. 
Most of the interest in this type of deformation in compact stars stems from the fact that it has 
the shortest growth time and leads to the emission of a strong gravitational-wave signal. 

The development of non-axisymmetric instabilities is commonly analyzed in terms of 
the quantity /3 = T/\W\ (i.e. the ratio between the kinetic rotational energy T and the 
gravitational potential energy W), that provides a dimensionless measure of the amount of 
angular momentum that can be tapped to feed the development of the instabilities. This 
parameter plays an important role in what is possibly the most celebrated of the non- 
axisymmetric instabilities: the so-called dynamical bar-mode instability. This is an m = 2 
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instability which takes place when the parameter (3 is larger than a critical one, j3 c . In the 
case of a Newtonian incompressible self-gravitating polytrope, for instance, the dynamical 
bar-mode instability develops for j3 > (3d = 0.2738 and is only weakly dependent on 
the considered polytropic index or whether the fluid is compressible. Post-Newtonian (PN) 
studies J4] or fully general-relativistic ones O correct this results only slightly, by reducing 
the threshold to somewhat lower values of the instability parameter. As an example, for a 
polytropic relativistic star with polytropic index r = 2, the accurate calculations reported 
in (6| reveal that the critical value is f3 c ~ 0.245 and that a simple dependence on the stellar 
compactness allows one to track this threshold from the Newtonian limit over to the fully 
relativistic one |7). 

The onset and development of the bar-mode instability has been traditionally studied by 
means of nonlinear 3D simulations of Newtonian stars that are either unmagnetized JS] [9] 
[TOj |TT] [T2l and, more recently, also magnetized lfl3l . In addition, PN and fully relativistic 
simulations have been performed and highlighted, for instance, that the persistence of the 
bar is strongly dependent on the degree of overcriticality and is generically of the order of 
the dynamical timescale. Furthermore, generic nonlinear mode-coupling effects between the 
m = 1 and the m = 2 mode appear during the development of the instability and these 
can severely limit the persistence of the bar deformation and eventually suppress the bar 
deformation j6j. These results have been recently confirmed by the perturbative calculations 
in Ql. 

Besides dynamical instabilities, which are purely hydrodynamical, secular instabilities 
are also possible in rotating compact stars and these are instead triggered by dissipative 
processes, such as viscosity or radiation emission. If, in particular, the dissipative mechanism 
is the emission of gravitational radiation, then the secular instability is also known as 
Chandrasekhar-Friedman-Schutz or CFS instability IfTSl [16l . Contrary to what their name 
may suggest, secular instabilities do not necessarily develop on secular timescales (although 
they normally do) and are characterized by having a much smaller threshold for the instability. 
Once again, in the case of a Newtonian polytrope, the critical secular instability parameter is 
as small as f3 c ~ 0.14 and thus much more easy to attain in astrophysical circumstances. 

Although widely observed in numerical simulations, the physical conditions leading to 
a dynamical bar-mode instability are difficult to be encountered in standard astrophysical 
scenarios. Such large values of the instability parameter, in fact, cannot be easily attained in 
old and cold neutron stars, which have been brought into uniform rotation and thus to rather 
small values of (3. However, more recently these pessimistic prospects have been changed 
when a new m = 2 instability has been discovered in differentially rotating Newtonian 
stars 0171 for values of /3 « 0.01, therefore well below the expected values for a dynamical 
bar-mode instability. The most salient aspect of this new instability is that it appears in stars 
with a large degree of differential rotation and that it grows on a timescale which is longer 
but comparable with the dynamical one. This instability has been referred to as the "low- 
TJ\ W\ instability" and its dependence on the polytropic index and on the degree of differential 
rotation has been studied in JT8). Since then, the instability has been observed or discussed 
in a number of related studies QU |20l ED |22l |23j |24] |25), all of which have highlighted the 
possible occurrence of this type of instability during the collapse of a massive stellar core. 

Despite the abundant numerical evidence on the development of this instability, the 
nature of these low-T/|V^| instabilities is still matter of debate and, most importantly, a 
sufficient criterion for its onset has not been derived yet. This instability has been studied in 
great detail by Watts and collaborators [26 , 2 1, who have made a number of phenomenological 
predictions either using a toy shell-model first introduced in 1271 1281 , or for a stellar model in 
Newtonian gravity. Overall, the work of Watts and collaborators (but see also l29l ) recognizes 
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the low-T/|W| instabilities as the manifestation of a more generic class of instabilities, the 
shear instabilities Q , that is unstable oscillations that do not exist in uniformly rotating 
systems and are associated to the existence of a corotation band ll30ll3Tl . Watts, Andersson 
and Jones suggest, in particular, that a necessary condition for the development of the 
instability is "corotation", that is the presence of a point at which the star rotates at the same 
pattern speed of the unstable mode ||2 |. An alternative suggestion on the necessary conditions 
has been made also by Ou and Tohline ll20l . who instead associate the development of the 
instability to the presence of a minimum in the vortensity profile of the star. This minimum 
can then drive unstable not only the corotating m = 2-modes but also the odd modes such as 
the m = 1 and m = 3-modes 11321 . In this interpretation, the growth time of the instability 
is proportional to distance between the corotation radius, i.e. the radial position at which the 
unstable mode corotates with the star, and the minimum of the vortensity. 

The purpose of this work is to shed some light on the development of shear instabilities 
and, in particular, to validate one prediction made, although not explicitly, by Watts, 
Andersson and Jones. More specifically, we show that, for sufficient amounts of differential 
rotation, shear instabilities develop for any value of the instability parameter j3 and also below 
the expected critical value for the dynamical bar-mode instability. We therefore provide 
evidence that the low-T/|VK| instability is not a new instability but rather the manifestation of 
a shear instability in a region where the latter is possible only for small values of /3. 

Our analysis proceeds via the simulation in full general relativity of sequences of neutron 
star models having constant rest-mass and constant degrees of differential rotation, but with 
different amounts of rotation, i.e. with different values of (3. The neutron-star matter is 
described by a realistic equation of state (EOS) defined by the unified SLy prescription HI 
and we study the development of the non-axisymmetric instabilities from their linear growth 
up to the fully nonlinear development and suppression. Interestingly, we find that depending 
on the degree of differential rotation, the shear instability leads either to the growth of a 
single modes (for the low-/? models) or to the simultaneous presence of up to three unstable 
modes (for the high-/? models), which produce beatings in the growth of the overall m = 2 
deformation. Special attention is also paid to the properties of the unstable modes and to 
their position within the corotation band or the vortensity profiles. In this way we are able 
to confirm both the necessary conditions proposed so far for the onset of the instability. In 
particular, we show that all the unstable modes are within the corotation band of the progenitor 
axisymmetric model (cf. j2 j) and that all of the unstable models have vortensity profiles with 
a local minimum (cf. Il20ll ). 

The structure of the paper is as follows: in section [2] we describe the numerical setting 
of our simulations, the EOS we used and the initial models we generated. In section [3] we 
describe the quantities and tools we used to monitor the evolution of the instability. In 
section [4] we report the results of the simulations and section [5] is dedicated to conclusions 
and discussion. We use a spacelike signature (—,+,+,+) and a system of units in which 
c = G = Mq = 1 (or in cgs units whenever more convenient). Greek indices are taken 
to run from to 3, Latin indices from 1 to 3 and we adopt the standard convention for the 
summation over repeated indices. 

2. Numerical Setup and Initial Models 

In what follows we provide a brief overview of the numerical setup used in the simulations, 
of the realistic EOS adopted and on the procedure followed for the construction of the initial 
axisymmetric models. 
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2.1. Numerical Setup 



We solve numerically the full set of Einstein equations 



(1) 



where G^ IU and T^ v are the Einstein tensor and the stress-energy tensor, respectively. 
The equations are solved within the "3+1" decomposition of spacetime, in which the 4- 
dimensional metric is decomposed into the spatial metric jij, the lapse function a and the 
shift vector components The field equations, which then also provide an evolution for the 
extrinsic curvature tensor Kij, are then coupled to those of general relativistic hydrodynamics 



Above w M is the fluid 4-velocity, p is the fluid pressure, e the specific internal energy and 
p the rest-mass density, so that e = p(l + e) is the energy density in the rest frame of the 
fluid. The set of hydrodynamics equation is then closed by a prescription for the properties 
of the matter in the form of a relation between the pressure and other quantities in the fluid, 
e.g. p = P(p, e), and for which we have chosen a cold and realistic EOS which will be 
discussed in the following section. 

The evolution of the spacetime was performed using the C CAT IE code, a three- 
dimensional finite-differencing code providing a solution of a conformal traceless formulation 
of the Einstein equations (see l33l for the explicit expressions of the equations solved in 
the code and also l34l for a more recent and improved implementation). The relativistic 
hydrodynamics equations, on the other hand, were solved using the Whisky code, which 
adopts a flux-conservative formulation of the equations as presented in [35] and high- 
resolution shock-capturing schemes or HRSC (see l36l [371 [38ll for the explicit expressions 
of the equations solved in the code and also ll39l for a more recent extension of the code 
to MHD). The Whisky code implements several reconstruction methods, such as Total- 
Variation-Diminishing (TVD) methods, Essentially-Non-Oscillatory (ENO) methods l40l and 
the Piecewise Parabolic Method (PPM) HP . Also, a variety of approximate Riemann solvers 
can be used, starting from the Harten-Lax-van Leer-Einfeldt (HLLE) solver l42l . over to 
the Roe solver J43j and the Marquina flux formula fl4l (see |[36l [371 for a more detailed 
discussion). All the results reported hereafter have been computed using the Marquina flux 
formula and a PPM reconstruction. 

Both the Einstein and the hydrodynamics equations are solved using the vertex-centered 
adaptive mesh-refinement (AMR) approach provided by the Carpet driver 1451 . Our rather 
basic form of AMR consists of box-in-box structures centered on the origin of the coordinate 
system and with the finest grid covering the whole star at all times. The simulations reported 
here make use of 4 levels of refinement, with the finest having a resolution of 221 m and the 
coarsest one a resolution of 1.77 km. The outer boundary was set relatively close to the star 
and at a distance of ~ 159.5 km, i.e. at about ~ 10 times the size of the star. A reflection 
symmetry across the (x, y) (equatorial symmetry) plane was used to reduce the computational 
costs, but not a rotational one around the z-axis (7r-symmetry) as it would have artificially 
prevented the growth of odd-m modes (see discussion in |6l). 




(2) 




(3) 
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2.2. Realistic Equation of State 

As mentioned above, the system of hydrodynamics equations needs to be closed by an EOS 
relating the pressure with the other primitive variables, e.g. the rest-mass density. Previous 
studies of the bar-mode instability, both in Newtonian gravity and in general relativity, have 
been focused on the use of ideal fluids and analytic EOSs, either in the form a of a "polytropic" 
(and isentropic) EOS p = p(p), or of an "ideal-fluid" and (non-isentropic) EOS p = p(p, e) 
(cf. discussion in 13). While these two descriptions are expected to provide results that 
are qualitatively correct, a more accurate modelling of these instabilities in compact stars 
necessarily requires a more physically-motivated description of the neutron-star matter. 

It is in this spirit that we have here considered a realistic EOS, namely the unified SLy 
EOS O, which models high-density and cold (i.e. zero temperature) matter via a Skyrme 
effective potential for the nucleon-nucleon interactions. The SLy EOS, which describes via 
a single effective Hamiltonian the neutron star's interior, is supplemented with the HP94 
EOS [46] to describe the crustal matter and with the BPS EOS BTll for lower density 
regions. This prescription results in a one-parameter EOS in the form p = p(c(p)) = p(p), 
where the SLy EOS is used for p > 4.979 x lO^gcm" 3 , the HP94 EOS is used for 
I0 8 gcm- 3 < p ^ 4.979 x 10 w gcm~ 3 and the BPS EOS for p < 10 8 gcm" 3 (see also 
figure 1 of [48 ]). In addition, at even lower densities the EOS becomes temperature dependent 
(and thus no longer a simple barotropic EOS), but because these these regions are well below 
the threshold for the artificial atmosphere, we do not consider an additional prescription for 
P iS 10 8 g cnr 3 . 

We recall, in fact, that our HRSC methods require the use of a tenuous atmosphere which 
fills the regions of the computational domain not occupied by the compact star. The threshold 
value for the rest-mass density of the atmosphere is chosen to be several orders of magnitude 
smaller than the maximum value and in our simulation a fluid element is considered to be part 
of the atmosphere if its rest-mass density p satisfies pj max (p(t = 0)) < x 10~ 8 . When this 
happens the fluid element is treated as a non-dynamical cold fluid described by a polytropic 
EOS, p(p) = Kp T , with r = 2 and K = 100 and its velocity is set to zero (see l49l for a 
more detailed discussion on the use of the atmosphere in the Whi sky code). 

The practical implementation of the realistic EOS can take place in a number of different 
ways. The simplest is to use standard interpolation techniques, e.g. based on Lagrangian 
polynomials, on the values of the published tables. While straightforward, the interpolations 
in this approach do not guarantee in general that the thermodynamics relations are fulfilled 
(see l50l for a thermodynamical preserving interpolation). In addition, the derivatives of 
the fields, e.g. of the pressure to evaluate the sound speed, are typically not available in the 
tables. Furthermore, the use of high-order interpolation and/or finite differences can lead to 
undesirable spurious oscillations. 

A second approach that removes all of these problems, uses analytic fits that have been 
proposed for the pressure. As an example, ref. 148) suggested to fit the specific internal energy 
of the unified SLy EOS table with the expressiory 

6 = %7 P l^ k{ - pMp)+Pl)} + 8xlO-a + 2.1p0.58 5 M^(log(p) +Pr)} ,(4) 

where fo{x} = l/(e x + 1), p and e are in cgs units, and the coefficients pi are p- t = 
{0.320, 2.17, 0.173, 3.01, 0.540, 0.847, 3.581} (see Table 2 of |48)). Equation © is 
obtained from Eq. (15) of |48"1 after using p = m^n, where n is the baryon number density 
and me = 1-66 x 10~ 24 g is the mass of the nucleons. As discussed in ifBTl . it is then possible 

| Note that a different notation is used in 1481 for some primitive variables. 
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to compute the pressure from the value of e using the first principle of thermodynamics at 

T = 

9 de 

P = P-r, (5) 
dp 

and thus to have an evaluation of the pressure which is thermodynamically consistent. 
The differences between the fit and the table are typically less then 2%. Unfortunately, 
although apparently very convenient, the evaluation of the fitting formulas containing several 
exponential and logarithmic functions, turns out to be computationally rather expensive even 
if done in a optimized way. 

As a third approach, which combines the efficiency of a table search with the 
thermodynamical consistency of an analytic fit, consists of performing a simple linear 
interpolation among the tabulated values constructed from the analytic fit. Besides being 
highly efficient, a linear interpolation also eliminates the spurious oscillations that arise, for 
instance, in the derivative of the pressure if high-order interpolation formulas are used. In this 
case, the interpolation error can be reduced simply by populating the analytically constructed 
tables with a large number of entries, e.g. ~ 600 in place of the ~ 150 which are typically 
available in published tables. This third approach is the one actually implemented in Whisky 
and provides a speed up of about 20% with respect to the evaluation of the pressure via the 
analytic fits and with comparable accuracy. 



2.3. Initial Data 

The initial data for our simulations are prepared as stationary and axisymmetric equilibrium 
solutions for rapidly rotating relativistic stars [52). Adopting spherical quasi-isotropic 
coordinates, the line element of the corresponding spacetime is 

ds 2 = -c» +u dt 2 + c»-"r 2 sin 2 9(d(f> - udt) 2 + c 2 «(dr 2 + r 2 d9 2 ) , (6) 

where p, v, u> and £ are functions of r and 9. Moreover we assume the usual relativistic 
j-constant law of differential rotation and that amounts to assume an angular-velocity 
distribution of the form 



n r - n 



A 2 r 2 



(CI - uj)r 2 sin 2 9e 



-2// 



1 - (fl - w) 2 r 2 sin 2 9e~ 



(7) 



where r e is the coordinate equatorial stellar radius and the coefficient A provides a measure of 
the degree of differential rotation. Expression (0 represents the general-relativistic equivalent 
of the simpler Newtonian j-constant law fl27l 

(A 2 r 2 + r 2 sin 2 9) 

Clearly, A ~ > oo corresponds to a star in uniform rotation, while A — >• corresponds to a star 
with increasing degree of differential rotation. As a reference, A = 1 yields a star with an 
angular-velocity profile which varies of a factor ~ 3 between the center and the surface of the 
star (cf. left panel of figure|2]). 

In practice, we have computed a very large number of initial models using the SLy 
prescription for the EOS for which we have computed baryonic mass Mb, the gravitational 
mass M, the angular momentum J, the rotational kinetic energy T and the gravitational 
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Figure 1. Left panel: position in the (/J, M/R e ) plane of the initial models computed with 
A = 1, with the filled circles representing those we have evolved numerically (A similar 
behaviour is shown also by the models with A = 2. Indicated with solid thin lines are 
isocontours of constant baryon mass models while indicated with a thick dashed line is the 
threshold to the dynamical bar-mode instability as computed for a T = 2 polytrope (7). Note 
that the threshold for the instability tends to increase for smaller rest-masses. Right panel: the 
same initial models as in the left panel but shown in {pc,r p /r e ) planes with isocontours of 
constant ft (upper part) or constant M/R e (lower part). See TablefJJfor a summary of the 
properties of the initial models. 



binding energy W defined as 

M b = J d 3 x ^W L p, M= Jd 3 x (-2T ° + T£) a ^7 , (9) 

E int = Jd 3 x^W L pe, J = Jd 3 xT°a^, (10) 

T=^Jd 3 xnT°a^, W = T + E int + M b -M, (11) 

where y^y is the square root of the determinant of three-dimensional metric 7^ and W L = au° 
is the fluid Lorentz factor. We stress that the definitions (l9l-(fTTTi of quantities such as J, T, W 
and (3 are meaningful only in the case of stationary axisymmetric configurations and should 
therefore be treated with care once the rotational symmetry is lost. 

Out of this large set, we have then selected for the numerical evolution a number of 
models so as to build sequences of constant baryonic mass with Mj = 2.5 M Q , degree of 
differential rotation A = 1,2 and values of the instability parameter f3 ranging between 

0. 140 and 0.250. Each model in this sequence is supramassive, namely it has a mass 
which is larger than the maximum mass allowed for a corresponding nonrotating model, 

1. e. M ma x\a=o = 2.05 M Q , M b%max \ n =o = 2.43 M Q , although it is not hypermassive, 
namely it does not have mass which is above the maximum mass for a uniformly rotating 
model, i.e. M ma , x \ n =n max = 2.41 M , M 6 , moa; |n=n mM = 2.84 Mq. There are two 
different reasons why such a large-mass model has been chosen. The first one is that we 
are interested in the development of a shear instability in the metastable star produced by the 
merger of a binary system of two neutron stars. As shown by a number of authors and most 
recently in (49), the product of this merger is either a supramassive or a hypermassive neutron 
star. Hence, our reference model has a mass which is sufficiently large so as to be a reasonable 
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Table 1. Main properties of the simulated stellar models. Starting from the left: the name of 
the model, the differential rotation parameter Athe instability parameter /3, the compactness 
M/R e , the central rest-mass density p c , the maximum of the density p ma x (Note that for 
models with A = 2, p ma x = pc-), the ratio between the polar and the equatorial coordinate 
radii r p /r e , the proper equatorial radius R e , the gravitational mass M, the total angular 
momentum J divided by the square of the gravitational mass. 





Model 


A 




M/Re 


Pc/10 15 


Pmax/10 15 


r p /r e 


Re 


M 


J/M 2 














(g/cm 3 ) 


(g/cm 3 ) 




(Km) 


(M Q ) 




M 


. 1 


.140 


1 


0.140 


0.245 


0.989 


1.011 


0.562 


13.20 


2.19 


0.681 


M 


. 1 


.150 


1 


0.150 


0.240 


0.941 


0.968 


0.536 


13.48 


2.20 


0.708 


M 


. 1 


.160 


1 


0.160 


0.237 


0.893 


0.926 


0.511 


13.81 


2.21 


0.735 


M 


. 1 


.180 


1 


0.180 


0.226 


0.802 


0.847 


0.461 


14.49 


2.22 


0.789 


M 


. 1 


.200 


1 


0.200 


0.217 


0.712 


0.773 


0.413 


15.27 


2.24 


0.844 


M 


. 1 . 


.210 


1 


0.210 


0.212 


0.668 


0.737 


0.389 


15.70 


2.25 


0.873 


M 


. 1 . 


.220 


1 


0.220 


0.207 


0.618 


0.703 


0.365 


16.16 


2.26 


0.902 


M 


. 1 . 


.230 


1 


0.230 


0.200 


0.576 


0.669 


0.340 


16.67 


2.26 


0.933 


M 


. 1 . 


.241 


1 


0.241 


0.194 


0.522 


0.633 


0.312 


17.29 


2.27 


0.968 


M 


. 1 . 


.244 


1 


0.244 


0.193 


0.506 


0.624 


0.305 


17.47 


2.28 


0.978 


M 


. 1 . 


.247 


1 


0.247 


0.190 


0.490 


0.614 


0.297 


17.66 


2.28 


0.988 


M 


. 1 


.250 


1 


0.250 


0.188 


0.474 


0.604 


0.289 


17.85 


2.28 


0.999 


M 


.2 


. 125 


2 


0.125 


0.241 


1.143 


1.143 


0.642 


13.22 


2.17 


0.668 


M 


.2 


.150 


2 


0.150 


0.227 


1.039 


1.039 


0.578 


14.15 


2.17 


0.739 


M 


.2 


. 175 


2 


0.175 


0.210 


0.947 


0.947 


0.512 


15.33 


2.18 


0.801 


M 


.2 


.200 


2 


0.200 


0.184 


0.865 


0.865 


0.435 


17.34 


2.15 


0.878 



approximation to the product of a binary neutron star mergei|§ The second reason is that a 
sufficiently massive model is necessary in order to reach values of the instability parameter 
which are above the expected threshold for dynamical bar-mode instabilities as computed 
in Q, i.e. /3 max ~ 0.25. Indeed, the maximum possible value for f3 within the computed 
sequence is around /3 max — 0.2533 and thus just above the threshold (cf. long-dashed line in 
the left panel of figure[TJ. Note also that the threshold for the instability tends to increase for 
smaller rest-masses. 

Overall, when comparing with equilibrium models generated with a polytropic EOS with 
K = 100 and r = 2 (see table 1 of (6)) the realistic EOS models reach higher compactness 
(models in [;6 1 typically have Mj R e ~ 0.1) but lower values of (3 {i.e. /3 max ~ 0.2533 for the 
models considered here, while /3 max ~ 0.28 for the polytropic models considered in (6)). 

The whole space of parameters is shown in figureQ] whose left panel reports the position 
in the (/3, Mj R e ) plane of the initial models computed with A = 1, and where the filled circles 
represent those we have evolved numerically. Indicated with solid thin lines are isocontours 
of constant baryon mass models while indicated with a thick dashed line is the threshold 
to the dynamical bar-mode instability as computed for a T = 2 polytrope J3 ■ The right 
panel reports the same initial models considered in the left one but shown in (p c ,r p /r e ) 
planes with isocontours of constant j3 (upper part) or constant M/R e (lower part). The main 
properties of the simulated models are also reported in table [TJ where we also introduce our 
naming convention. Any initial model is indicated as M. % . #, with % being replaced by the 

§ We note that the use of a supramassive model has also the inconvenient consequence that it is not possible to 
construct a corresponding nonrotating model and this prevents us, for instance, from computing the frequency of the 
fundamental mode and compare it with the results of perturbative analyses. 
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value of the differential-rotation parameter A and # by the instability parameter (5. As an 
example, M . 1 . 2 is the star with A = 1 and f3 = 0.200. 

Finally, shown in figure [2] are the angular- velocity profiles (left panel) and the rest-mass 
density profiles (right panel) of some representative models, namely M. 1.150, M. 1.200 
and M . 1 . 250. Indicated with different symbols, which match the ones in the left panel 
of figure [5] are the normalized radial positions of the corotation radii, with the one for 
model M. 1 . 150 being shown filled to help distinguish it from the others. Although we 
postpone to sections |4~T1 and I4.2l the discussion of the implications of these corotation radii, 
two aspects of the initial data are worth emphasizing. The first one is that the amount of 
differential rotation for a given value of A effectively decreases when increasing the instability 
parameter j3 (cf. left panel of figure |2), thus resulting in a smaller corotation band for models 
with large j3. The second one is that all the initial models evolved are axisymmetric but have 
a "toroidal-topology", namely have the maximum density p max that is not at the center of the 
star, and thus p c < p max . This toroidal deformation increases with the rotation and thus with 
(i. 




o 0.5 




Figure 2. Left panel: Initial angular-velocity profiles for three representative models with 
small, medium and high values of and A = 1 (A similar behaviour is shown also by the 
models with A = 2). Indicated with different symbols, which match the ones in the left panel 
of figure [3] and [6] are the normalized radial positions of the corotation radii, with the one for 
model M . 1 . 1 5 being shown filled to help distinguish it. Right panel: the same as in the left 
panel but for the initial rest-mass density. 



3. Methodology of the analysis 

A number of different quantities are calculated during the evolution to monitor the dynamics 
of the instability. Among them is the quadrupole moment of the matter distribution 

I jk = J d 3 x ^W L px 3 x k . (12) 

which we compute in terms of the conserved density y^y W L p rather than of the rest-mass 
density p or of the Too component of the stress energy momentum tensor. Of course, the 
use of y^y W L p in place of p or of Too is arbitrary and all the three expressions would have 
the same Newtonian limit, though with different amplitudes for the gravitational waveforms 
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produced (see for example ll53lD . However, we here adopt the form (TT2l i because y/jW L p 
is a quantity whose conservation is guaranteed by the form chosen for the hydrodynamics 
equations. 

The quadrupole moment (TT~2T > can be conveniently used to quantify both the growth 
time of the m = 2 instability T2 and the oscillation frequency once the instability is fully 
developed a\. (Hereafter we will indicate respectively with 73 and af n \ the growth time and 
frequencies of the m = n unstable modes and we note that, as will be discussed later on, 
during the simulation a number of different modes appear, thus justifying the use of the upper 
index "i"). In practice, we perform a nonlinear least-square fit of the xy component of the 
computed quadrupole P (£) and we generally use as fitting function a sum of N (usually 
three) exponentially modulated cosines 

N 

P k (t) = ]T e*/™ cos(27r a\ n) t + , (13) 

i=l 

where P Q k = I^ k {t = 0). Because we commonly have only about 10 cycles in the time 
interval considered for the fits, extreme care needs to be applied when computing the growth 
time, especially when the oscillation frequencies and the growth times are close to each other. 
In these cases, in fact, variations of the initial phase of the modes <f>u-s can result in large 
variation of the growth times. In view of this, we will not report them. 

Using three components of the quadrupole moment in the [x, y) plane we can define the 
distortion parameters 77+ (t) and 77 x (t), as well as the axisymmetric mode r)o(t) as 

P x {t) - I*v(t) _ 2py{t) _ P x (t) + iw(t) 

v+ ^> = i*x(o) + iw(o) ' 1]x ^ = p x (o)+ivv(o) ' = P x (o) + iyy{o) ' 

so that the modulus t](t) and the instantaneous orientation of the bar are given by 

V (t) = W + (t) 2 +V,(t) 2 , = tan- 1 [ Ixx ^Z { % {t) ) ■ (15) 

Finally, as a useful tool to describe the nonlinear properties of the development and 
saturation of the instability, the rest-mass density is decomposed into its Fourier modes P m (t) 
as 

P m {t) = J<Pxpe im<1, , where </> = tan" 1 ^). (16) 

The phase <p m = arg(P m ) essentially provides the instantaneous orientation of the ?n-th mode 
when the corresponding mode has a nonzero power. Note that despite their denomination, the 
Fourier modes ( TTol l do not represent proper eigenmodes of oscillation of the star. While, in 
fact, the latter are well defined only within a perturbative regime, the former simply represent 
a tool to quantify, within the fully nonlinear regime, what are the main components of the 
rest-mass distribution. As a final comment we note that while all quantities (TT2T >— (fl~6b are 
expressed in terms of the coordinate time t and are not invariant measurements, the lengthscale 
of variation of the lapse function at any given time is always larger than twice the stellar radius 
at that time, ensuring that events on the same timeslice are also close in proper time. 



4. Results 



In what follows we first describe the dynamics of the shear instability as deduced from the 
numerical simulations and then contrast our results with the phenomenological predictions on 
the necessary conditions for its development. 
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Figure 3. Summary of the dynamics of some representative models, i.e. M. 1.150, 
M. 1.160, M. 1.200, M. 1.220, M . 1 . 2 4 1 and M . 1 . 2 5 , with increasing values of 
0, For each panel the upper part reports the evolution of the bar-distortion parameter 77+, 
while the lower part shows the evolution of the power in the different modes of the Fourier 
decomposition of the rest-mass density P m . See text for details. Note that we here report only 
some models and with A = 1 as they are representative also of those with A = 2. 
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Figure 4. Power spectral density (filtered using Harming windowing) in arbitrary units of 
the evolution of r)+ for the models described in figure [3] (black continuous line) and of the 
axisymmetric mode r)Q (red dashed line). Indicated with vertical lines (long-dashed black, 
dotted blue and dot-dashed magenta) are the peak frequencies reported in the left panel of 
figure|3]within the corotation band and which are used in figure|6]to mark the position of the 
corotation radii. Finally, shown with a shaded rectangular area is the corotation band for all 
models. 

4.1. Dynamics of the instability 



As mentioned when discussing the initial data, we have evolved numerically two sequences 
of constant baryonic mass, with M& = 2.5 Mq. The first sequence has a higher degree 
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of differential rotation {i.e. A = 1) and instability parameter ranging from j3 = 0.140 to 
/3 = 0.250. The second sequence, on the other hand, has a smaller degree of differential 
rotation (i.e. A = 2) and instability parameter between (3 = 0.125 and j3 = 0.200. For both 
sequences the highest value of f3 considered is also very close to the highest attainable with 
our initial data code (see figure [T). To keep the computational costs to an affordable level 
we have evolved of all these models for ~ 15 ms and although this time window is in general 
insufficient to capture the suppression of the instability, is adequate to measure the frequencies 
of the unstable modes and provide a first estimate of the growth times. 

After analyzing the results following the method outlined in the previous section and 
focusing on the properties of the distortion parameters 77+, t] x and 770, we find that all of the 
models show anm = 2 instability. The maximum distortions obtained during the simulation 
time are 77 ~ 0.1 (i.e. distortions of 10% with respect to the axisymmetric progenitor), with 
the maximum values being reached for the models of the A = 1 sequence. The models of the 
A = 2 sequence, in fact, have in general much smaller distortions, with rj ~ 0.01 over the 
timescale over which they have been evolved. We believe this is simply the consequence of 
the fact that the models in this sequence have smaller growth rates (see also the discussion in 
the next section). 

To show the general behaviour of the simulation we focus on the description of the 
result of the simulation for the A = 1 sequence as this is representative also of the one 
with A = 2. In particular, in the upper part of each panel in figure [3] we show the 
time evolution of the distortion parameter r/ + computed with Eq. ( TPfl i for six representative 
models: M.1.15 0, M.l. 160, M. 1.200, M. 1.22 0, M. 1.241 and M.l. 250. Similarly, 
in the lower parts of each panel we show the corresponding evolution of the Fourier-modes 
computed from Eq. ( TToT l for m = 1, 2, 3, 4. 

Note that for models with lower values of f3 the bar-mode deformation is very similar to 
the one already discussed in (6][7), growing exponentially and with only one unstable mode 
appearing (cf. first row of figure 0. However, as the rotation is increased, the development 
of the instability is more complex and at least two unstable modes appear which develop 
in different parts of the star. These two modes have very similar but distinct frequencies 
and growth times, leading to a series of beatings in the evolutions of whose irregular 
evolution makes the calculation of the growth times challenging (cf. second row). As the 
rotation is increased towards the maximal values of j3, up to three distinct modes appear and 
the evolution of the instability is correspondingly more complex (cf. third row). A similar 
behaviour is shown also by the evolution of the Fourier modes. The two models with lower j3, 
in fact, show a clear growth of the m = 2 mode and of the m = 1, with the latter becoming 
first comparable and then larger when the instability is suppressed. In models M . 1 . 2 5 , 
M. 1 . 241, on the other hand, the m = 2 and m = 3 modes have comparable amplitude 
for a long period and then the m — 2 becomes the dominant one (this is more clear in 
model M. 1 . 241). Finally, for larger values of (3, the m = 1 mode never attains values 
comparable with either the m = 2 or the m = 3, which instead control the evolution. 

Figure |4] reports the power spectral density (PSD) of r/ + and 770 for the same 6 models 
(solid black lines) and allows to appreciate how the spectrum changes as the instability 
parameter is increased. More specifically, it is very apparent that at low (3 the spectra have 
only one peak, whose maximum is marked by the vertical blue dotted lines and that is present 
at all the values of beta considered; this is the mode we refer to as erf (cf. table [2] and 
figure [3). As f3 increases, the amplitude of this peak decreases and it becomes the weakest 
for very large rotation rates. Starting from model M. 1 . 2 00 (although a hint of a peak is 
present already in M . 1 . 1 60), a second peak appears at higher frequency and is marked with 
magenta vertical dot-dashed lines; this is the mode we refer to as a\ and although it never 
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Table 2. Results of the analysis of quadrupole evolutions. The frequencies <r . and a\ are 
obtained from the position of the peaks in the PSD while the frequencies cr| are obtained 
from the nonlinear fitting of Eq. 1 1 3i . Reported in the last two columns are the edges of the 
corotation band as expressed in terms of angular frequency at the surface £l e /7r = Q(r = 
R e )/n and on the axis Q. c /tt. Note that depending on the rate of rotation and the degree of 
differential rotation some frequencies may not be present. 
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becomes the largest one, its amplitude increases with (3. Finally, a third peak appears at low 
frequency in the spectra starting from model M . 1 . 2 2 (even though a hint is present also for 
model M . 1 . 2 1 0) and is marked with black vertical long-dashed lines. This mode, which we 
refer to as a\, becomes the dominant one at high j3. Note also that in each panel of figure [4] 
we report with a dashed blue box the corotation band that will be further discussed in the 
following section. 

As a final remark we note that the peaks in the spectrum of I xx + I vy (red dashed 
lines in figure [U are related to axisymmetric m = modes and most likely to /-modes 
oscillation excited by the development of the instability. Indeed, their frequencies a f 
match reasonably well the phenomenological fit for the /-mode frequencies of nonrotating 
neutron stars with realistic EOSs computed in lf54l . Because our sequences do not contain 
nonrotating configurations (the rest mass is larger than the maximum one), this association is 
just qualitative and a more detailed investigation of the frequency spectra of the equilibrium 
models is necessary to confirm this suggestion. 



4.2. Necessary conditions for the instability 

To support the interpretation of these instabilities as shear instabilities we have considered 
whether the necessary conditions suggested by |2] and by l20l are met. We recall that Watts 
and collaborators pointed out that an unstable configuration should have the unstable mode 
with a frequency within the corotation band. In the case of the Newtonian expression ([8j, the 
corotation band is simpler to compute and for a mode with frequency a and azimuthal number 
m, this is simply given by 

VcA 2 <r _ 

< — <fl c . (17) 

A 2 + l m 
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We have therefore checked whether any of the unstable m = 2 (bar-modes) with computed 
frequency er| 2 N has pattern speed velocity, a\/m = cr|/2, within the corotation band. 
A careful and rather involved analysis has indeed confirmed the prediction of Watts and 
collaborators: namely, all the unstable modes are within the corotation band of the progenitor 
axisymmetric model. These results for the models with A = 1 are listed in table |2j where 
we report the edges of the general-relativistic corotation band in terms of angular frequency 
at the surface f2 e / 7r = £l(r = R c )/tt and on the axis il c /ir, as well as the frequencies of the 
unstable m = 2 modes obtained from the position of the peaks in the PSD (i.e. a l 2 ) or from 
the nonlinear fitting of Eq. ( TT3l l (i.e. a\). Note that depending on the rate of rotation and the 
degree of differential rotation some frequencies may not be present in the corotation band. 

The data in tablets also shown in the left panel of figure[5] where we plot the position of 
the non-axisymmetric frequencies a\ within the corotation band as a function of /?; indicated 
with a filled symbol, to distinguish it from the others, is the o\ frequency for model M . 1 . 1 5 
(cf. figure [2] and [3J. The right panel of the same figure shows the corresponding information 
and on the same scale but for the models with A = 2 (the inset shows instead a magnified 
view). Clearly, for both sequences all the unstable modes are within the band and, as the 
stellar rotation rate increases, more unstable modes appear for the same value of f3. 



1 

n /tt 
















p '////.■ 
- a~ ; 








- D ■'" □ 


□ 


:2 v D o 

O g 


CCCD 


I 




A 

i , , , 


= 1 ; 


0.15 




0.2 


0.25 



X 4 



n /tt 



n „/tt 




0.15 0.2 



A = 2 



ft 



Figure 5. Left panel: Position of the non-axisymmetric frequencies for models with A = 1 
and shown as a function of 0; indicated with thick solid lines are the edges of the corotation 
band, i.e. the frequency interval between f2 c /7T and Q e /ir. Shown with a filled symbol, to help 
distinguish it from the others, is the a\ frequency for model M. 1 . 150 (cf. figuref2]and|6j. 
Right panel: the same as in the left one and with the same scale but for models with A = 2; 
the inset shows instead a magnified view. Clearly, in both panels all the unstable modes are 
within the band and up to three unstable modes appear for the A = 1 sequence with increasing 



Watts, Andersson and Jones also give a qualitative description, summarized in their figure 
2 of 0, of how the unstable and stable models should be distributed in the (A, /3) plane. The 
considerations they make are particularly simple: for high degrees of differential rotation, that 
is at low values of A, the corotation band is rather wide and there will be both an upper value 
and a lower value of j3 between which the shear instability can develop (these critical values 
of j3 correspond to the entrance and exit of the unstable mode in the corotation band). This 
situation corresponds therefore to the one commonly encountered in numerical simulations, 
such as those in |[T8l[l^l2^l2~Tll22ll23ll24l . and for which the shear instability takes place only 
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for very small values of the instability parameter and on timescales that are much longer than 
the dynamical one. When moving to larger degrees of differential rotation, that is when going 
to smaller values of A, the corotation band becomes larger and larger and the shear instability 
can develop essentially for all values of (3, merging with the dynamical bar-mode instability 
for f3 > 0.25. This is exactly what has been found here. Conversely, when moving to smaller 
degrees of differential rotation, that is when going to higher values of A, the corotation band 
becomes thinner and the shear instability can develop only for a smaller range of f3. As the 
differential rotation is further decreased and the star tends to rotate uniformly, the corotation 
band width vanishes, all the models are stable to the shear instability and subject only to the 
dynamical bar-mode instability for (3 > 0.25. This is indeed the case for the unstable models 
evolved in |6|, which were purely (bar- mode) dynamically unstable and none of which had 
the unstable mode within the corresponding corotation band, but above it. 

In essence, therefore, there should be an intermediate range of A for which the instability 
is absent at low f3, appears at intermediate values and then disappears again at high /?, thus 
defining an interval of values of f3 for which the models are unstable. To validate also this 
prediction we simulated a sequence of models with smaller differential rotation and A = 2. 
This seconds sequence has the same baryonic mass as that with A = 1, but with (3 in a smaller 
range, namely between j3 = 0.125 and j3 = 0.200 (Note that j3 = 0.200 is also very close 
to the largest value for which we could build an equilibrium model.) Unfortunately, also all 
of these models show an m = 2 shear instability, with the unstable frequencies falling within 
the corotation band, as it was for the A = 1 stars. Of course, lack of evidence is not evidence 
of absence and the fact that we have not found stable models within our range of values of f3 
most likely means that these stable models have to be searched either for values of j3 < 0.125 
(for this value of A = 2), or by moving towards higher values of A, where the corotation band 
is less large. 

An obvious consequence of the phenomenological scenario described by Q and 
confirmed by the calculations reported here is that the idea of a low-T/|VK| instability 
is indeed misleading and it is instead more meaningful to think of a more generic shear 
instability that, depending on the degree of rotation and of differential rotation, may manifest 
itself on timescales that are comparable with the dynamical ones (as in the cases reported 
here) or on much longer ones (as in the cases reported in ifTTll ). 

Another phenomenological prediction made by Watts, Andersson and Jones in (2 ||, is that 
the growth times should be shorter in the center of the band and increase towards the edges (the 
growth times should in fact diverge at the edges). As shown in figure |5j all of the simulated 
models do have unstable modes well inside the corotation band and this probably explains 
why we observe them develop on a timescale which is comparable with the dynamical one 
and not on much longer timescales as in the original finding of (T71 . We also find indications 
that the A = 2 models, which generically have longer growth times, have unstable modes 
which occupy regions of the corotation band that are overall more central (cf. right panel of 
figure|5]l and thus in contrast with what expected from 0. However, the difficulties mentioned 
above in computing an accurate estimate of the growth times and the intrinsic difficulties in 
determining what is the central part of the corotation band, prevent us from providing a more 
quantitative validation of this prediction. 

We now switch to consider our results in the light of the other necessary condition for 
the onset of the shear instability discussed in l20l . As mentioned in the Introduction, it has in 
fact been shown that the presence of a minimum in the profile of the vortensity is a necessary 
condition for a mode in corotation to be unstable l20l . The intuitive description is that the 
vortensity well can act as a resonant cavity inside the star, amplifying the modes that happen 
to lay near its minimum l32l . Indeed, the growth rate of the unstable mode is expected to 
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Figure 6. Left panel: radial profiles of the vortensity V for three representative models with 
small, medium and high values of f). The different symbols match the ones in figure[2]and in 
the left panel of figure [5] with the one for model M. 1 . 150 reported as filled, and show the 
actual position of the corotation radii. Note that none of these coincides with the minimum of 
V; for compactness we have reported only the models with A = 1. Right panel: normalized 
corotation radii for the different frequencies a\ presented in the left panel of figure[5]shown as 
a function of /3; the same convention of the previous figures is used for the different symbols. 



depend on the location of its corotation radius with respect to the vortensity profile, being 
proportional to the depth of its corotation radius inside the vortensity well [20|. Note that the 
mere existence of a local minimum in the vortensity cannot be used as a sufficient condition 
for the occurrence of a shear instability. All of the unstable models considered in (6), in 
fact, do show a local minimum in the vortensity but do not have the unstable modes in the 
corotation band. 

To validate whether this Newtonian condition holds also for the general-relativistic 
instabilities simulated here, we have computed for all the models the Newtonian vortensity, 
which is defined as the ratio, along the radial cylindrical coordinate, between the radial 
vorticity and the density 

V=^, (18) 
P 

where w is the radial cylindrical coordinate (a fully general-relativistic definition of the 
vortensity is also possible but more complicated to compute and not significantly different 
from the Newtonian one). When doing this we found that all the unstable models have 
vortensities with a local minimum in the star and this is shown in the left panel of figure 
[6] which reports the radial profiles of the vortensity V for three representative models with 
small, medium and high values of /3. The different symbols match the ones in the left panel 
of figure [2] and [5] with the one for model M . 1 . 1 5 reported as filled, and show the actual 
positions of the corotation radii. The right panel of the same figure, on the other hand, shows 
the normalized corotation radii for the different frequencies er| presented in the left panel of 
figure[5]as a function of /3. 

Note that none of these coincides systematically with the minimum of V, nor with the 
maximum of the rest-mass density. All of the corotation radii, however, do move towards 
larger radial positions as the rotation rate is increased, exactly as does the minimum of the 
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vortensity and the rest-mass maximum. Although we cannot confirm that modes whose 
corotation radius is closer to the minimum of the vortensity have systematically shorter 
growth times (our data is not sufficiently accurate for this), we can compare the models 
with the higher degree of differential rotation A = 1 with those having a smaller degree 
of differential rotation A = 2. In the first case the vortensity well is considerably deeper, with 
AV/V = 1 — Vmi n /V r= o ~ 0.6 — 0.7, while for the latter the vortensity has shallower wells 
with AV/V = 1 - V min /V r=0 ~ 0.02 - 0.25. Since the models with A = 2 have smaller 
growth rates, our results indicate therefore that also the existence of a local minimum in the 
vortensity can be taken as a necessary condition for the development of the shear instability 
and that the depth of the vortensity well can be used to estimate the growth of the instability. 

5. Conclusions 

For many years the properties of rapidly rotating and self-gravitating fluids have been 
characterized by a complete analytic perturbative theory which provided, for instance, 
sufficient conditions for the development of instabilities (see, for instance, [3] [55] for a 
collection of results). As numerical simulations have become increasingly accurate and 
stable on longer timescales, these predictions, both in Newtonian theory and in full general 
relativity, have been verified, corrected and in some cases extended. Our ability of modelling 
such configurations has now reached a maturity such that a number of new properties and 
instabilities have been "discovered" numerically, but which do not have behind a fully 
perturbative description. The main reason for this is that most of this phenomenology is 
the result of physical scenarios which are much more complex than the ones investigated 
perturbatively in the past, e.g. non-isolated systems with high differential rotation and 
exchanging mass and angular momentum, which are much more difficult to treat analytically. 

A most notable example of these complex and yet ubiquitous instabilities is the so- 
called "low-T/|W| instability", which was initially found in IfTTl and then reproduced in a 
number of other different scenarios j[T8] [19J [20j [21] [22l [23l [24] (25] . As the phenomenological 
description of this instability provided by the simulations has become richer and richer, a 
full understanding of the mechanisms that lead to its development has lagged behind and 
it is presently unclear. We are therefore in a situation in which numerical simulations can 
probe regimes and conditions which are not yet accessible to perturbative calculations, and 
can guide the latter by confirming or refuting those scenarios that although possible do not 
find a realization in practice. 

This work has followed this spirit and has used fully general-relativistic calculations of 
rapidly and differentially rotating neutron stars modeled with a realistic EOS to shed some 
light on the development of the low-T/|W| instability. In particular we have concentrated 
our attention on validating an indirect prediction, made by Watts, Andersson and Jones J2], 
who recognized the low-T/|W| instability as the manifestation of a more generic class of 
instabilities associated to the existence of a corotation band ll30l[3T1 . the shear instabilities, 
and that should develop for any value of the instability parameter f3 when sufficient amounts of 
differential rotation are present. This is exactly what we have found in our simulations. More 
specifically, we have performed simulations of sequences of neutron-star models described by 
a realistic SLy EOS (T) and having constant rest-mass and degrees of differential rotation, but 
with different amounts of rotation. In all cases considered we have found the development of 
a bar-mode instability growing on a dynamical timescale, even when the initial axisymmetric 
models were well below the critical limit for the dynamical bar-mode instability. These 
results, which match well the phenomenological scenario portrayed in (2), suggest therefore 
that the idea of a low-T/jW^I instability is indeed misleading and should be replaced by 
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the more general one of shear instability. Depending then on the degree of rotation and 
of differential rotation, the instability will develop on timescales that are comparable to the 
dynamical one (as reported here) or on much longer ones (as reported in the first low-T/|W] 
instability studies). 

Special attention has also been paid to the properties of the unstable modes and to their 
position within the corotation band or the vortensity profiles. In particular, we have shown that 
all the unstable modes are within the corotation band of the progenitor axisymmetric model 
(which is the necessary condition for the development of the instability proposed by J2)) 
and that all of the unstable models have vortensity profiles with a local minimum (which is 
the necessary condition suggested by ll20l ). Finally, by comparing the growth times among 
models with different degree of differential rotation we have shown that there is a correlation, 
although not a strong one, between the depth of the vortensity well and growth rate of the 
instability, with the latter being larger for models with deeper wells. 

In summary, the results presented here shed some light on several aspects of shear 
instabilities, but they also reveal that more work is required, for instance, to distinguish 
between the predictions based on the corotation band and the ones based on the vortensity 
well, or to establish whether in effect they just represent two different ways of expressing 
the same physical conditions. Clarifying these aspects will requires additional analytical and 
numerical modelling, and this will be part of our future research. 
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